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The general relativistic (Mercury-type) periastron advance is calculated here for the first time with exquisite 
precision in full general relativity. We use accurate numerical relativity simulations of spinless black hole bi- 
naries with mass ratios 1/8 ^ m\/m2 S5 1 and compare with the predictions of several analytic approximation 
schemes. We find the effective-one-body model to be remarkably accurate, and, surprisingly, so also the predic- 
tions of self-force theory [replacing rai/m2^>»2im2/(rai+»i2)^]- Our results can inform a universal analytic 
model of the two-body dynamics, crucial for ongoing and future gravitational-wave searches. 
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Introduction. The anomalous rate of Mercury's perihelion 
advance was originally recognized in 1859 by the astronomer 
Urbain Le Verrier. For the first time, Newton's law of uni- 
versal gravitation could not be reconciled with observation. 
Treating Mercury as a test body in free fall in the gravitational 
field generated by the mass Mq of the Sun, Einstein derived 
the lowest order (weak-field) general relativistic angular ad- 
vance per orbit [IJ 
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where a and e are the semi-major axis and eccentricity of Mer- 
cury's orbit, respectively. Equation ([TJ perfectly accounted for 
the observed discrepancy of ^ 43" per century, thus providing 
the first successful test of general relativity. More recently, the 
same effect — but with a much larger amplitude, of order a few 
degrees per year — has been observed in the orbital motion of 
binary pulsars |2 |. Today, the exciting prospects of observing 
gravitational waves from the inspiral and merger of compact 
binaries, using interferometric detectors like LIGO or Virgo, 
provide a modern context for the problem of relativistic peri- 
astron advance, and a motivation to go far beyond Einstein's 
weak-field test-particle approximation. 

In this Letter we restrict our attention to binaries composed 
of two black holes. Their orbital dynamics can be analyzed us- 
ing several approximation schemes in general relativity: post- 
Newtonian expansions |3|, black hole perturbation theory 
and the effective-one-body model [|5|. It can also be studied 
using fully nonlinear numerical relativity (NR). While NR can 
now routinely perform accurate binary black hole simulations 
0, approximation methods remain valuable given the high 
computational cost of these simulations, and their restricted 
utility when the mass ratio is too extreme. It is important to 
assess the predictions of the various approximations against 
the NR benchmark, since (i) it allows crucial cross-validation 
tests, (ii) it helps delineate the respective domains of validity 
of each method, and (iii) it can inform the development of a 
universal semi-analytical model of the binary dynamics. 



Neglecting radiation reaction, the motion of two non- 
spinning black holes on a generic eccentric orbit involves two 
frequencies: the radial frequency (or mean motion) Q.r, and 
the averaged angular frequency £2(p, defined by 
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where P is the radial period, i.e. the time interval between two 
successive periastron passages, <p — d<p/dt '\s the time deriva- 
tive of the orbital phase ^(f), and A«l>/(27r) = A' - 1 is the 
fractional advance of the periastron per radial period. In the 
circular orbit limit, the relation between K = Q.,p/Q.r and Q.ip 
is coordinate invariant (for a large class of physically reason- 
able coordinate systems), and therefore provides a natural ref- 
erence for comparing between the predictions of the analytical 
and numerical methods currently available. 

In this Letter we present new accurate NR simulations start- 
ing at lower orbital frequencies than in previous work L7 -9). 
We outline the respective computations of the invariant re- 
lation K{Q.(p) in numerical relativity, post-Newtonian theory, 
the effective-one-body formalism, and black hole perturbation 
theory. We then perform an extensive comparison which, for 
the first time, (i) encompasses all of these methods, and (ii) 
focuses on the orbital dynamics of the binary, rather than the 
asymptotic gravitational waveform. We also discuss the im- 
plications for the modelling of coalescing compact binaries. 
(We henceforth set G = c = 1.) 

Numerical Relativity. The periastron advance of non- 
spinning black hole binaries was estimated for the first time 
in general relativistic numerical simulations in fTOl. In the 
present work, we improve considerably on the accuracy of 
these calculations. Our results are based on new and longer 
simulations of the late stage of the inspiral of black hole bi- 
naries, using the Spectral Einstein Code SpEC lflTl[T2l . with 
mass ratios q = m\/m2 between 1:1 and 1:8, and eccentrici- 
ties e in the range [0.0015,0.023]. These runs are summarized 
in Table [l] and will be described in detail elsewhere ll9l [T3l . 
(Ref. lilOi discusses the definition of e in these simulations.) 
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1/3 
1/5 
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1/8 



19 0.021 34 0.9949 0.589 

18 0.023 27 0.9950 0.573 

14 0.002 29 0.9821 1.692 

14 0.008 23 0.9879 1.154 

13 0.015 20 0.9890 1.071 



-79.1 0.0111 0.0312 

-75.9 0.0129 0.0316 

-87.1 0.0181 0.0313 

-62.8 0.0183 0.0361 

-57.0 0.0204 0.0333 



13 0.0015 24 1.0028 -0.099 -26.8 0.0197 0.0355 



TABLE I: Simulation parameters. Here q = ini/ 1112, m = mi +1112, d 
is the initial coordinate separation, e the initial eccentricity, and A'oib 
the total number of orbits in the simulation. The fitting parameters 
{ao,"li"2} [cf- Eq. ([3|] are computed for the restricted frequency 
range fl, ^ f2,p ^ Q.f. 



We compute Q.(p and Q.r using the orbital frequency Q.{t) 
extracted from the motion of the apparent-horizon centers (in 
harmonic coordinates): let c, (f) be the coordinates of the cen- 
ter of each black hole, and define their relative separation 
r = Ci — C2; then f2 = |r x r|/r^, where the Euclidean cross 
product and norm are used. The frequency Q.{t) can be writ- 
ten as the sum of a secular piece (given by the average fre- 
quency Q.ip) and a small oscillatory remainder — ^both of which 
drift slowly in time due to radiation reaction. To compute 
/Tnr at some coordinate time T, we choose a time interval of 
width W X 27z/Q.{T), centered on T, and fit i2(f) to the model 
^{t) = PoiPi - tf^ + Pi cos [iH + ps{t-T)+ pe{t - r)2] , 
where the /7,'s are fitting parameters. We then write £l(f,{T) = 
Poipi — TY' and D.r{T) — ps, compute the ratio K}^r{T) = 
Clip{T) /Q.r{T), and hence obtain /Tnr as a function of Q.(p. Fi- 
nally, we fit /rNR(£2(p) to a smooth quadratic polynomial using 

■'^NR = [ao + ai{mD.^)+a2{mQ.^)^]Ks, 
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where m = mi +m2 is the total mass of the binary. The results 
of the fits are given in Table [l] For convenience, the numeri- 
cal periastron advance A^nr is normalized by the test-particle 
result A'schw, which is known in closed form as llT4l ITSi 
Kschw = (1 ^ 6x)^'/2, where x = (mil^)^/-' is the usual di- 
mensionless coordinate invariant post-Newtonian parameter. 

The variance in the numerical data for various window sizes 
W provides an estimate of the error in Kj^r. We point out 
that the finite (non-zero) eccentricity in the NR simulations 
introduces a small error, since we are interested in the e ^ 
limit. However, as the leading-order result ([TJ suggests, and 
calculations at higher post-Newtonian (PN) orders confirm, 
this error scales like e^, which in our simulations is always 
< 5 X 10^^, and decreasing monotonically with time. 

The numerical data form the basis for our comparisons. We 
will now discuss the different approximation schemes in turn, 
summarizing the results in Figs.[T]and|2](showing Kasa func- 
tion of frequency for two fixed mass ratios), and Fig. |3](show- 
ing K asa function of mass ratio for a given frequency). 

Post-Newtonian Theory. Einstein's result ([T]) was general- 
ized to arbitrary masses nii and m2 by Robertson |fT6ll . Fol- 
lowing the discovery of binary pulsars in the 1970s, an im- 
proved modelling of the orbital dynamics of these compact 



binaries was required, leading to the extension of this IPN re- 
sult to 2PN order 1 14 1. [As usual we refer to nPN as the order 
equivalent to terms ^(c^^") in the equations of motion be- 
yond the Newtonian acceleration.] More recently, the need for 
extremely accurate gravitational-wave templates modelling 
the inspiralling phase of coalescing compact binaries moti- 
vated the computation of the equations of motion through 3PN 
order. These results allowed also the calculation of the perias- 
tron advance at the 3PN accuracy for eccentric orbits ifTTl . 

For quasi-circular orbits, combining Eqs. (5.8) and (5.25) 
of Ref. 1 17 1, we obtain the 3PN-accurate expression of K as 
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+ lv^]x^ + i)'{x'^). (4) 



The symmetric mass ratio v = mini2/ni^ is such that v = 1/4 
for an equal mass binary, and v — > in the extreme mass ratio 
limit. The term °^ in Eq. Q, which is a 3PN effect, con- 
tributes less than 1 % to Kjpf^, for all mass ratios. This suggests 
that the exact K may be well approximated by a linear func- 
tion of V. Figures [TJ|3] show a good agreement between the 
3PN and NR results for g = 1, with < 1% relative difference 
even at the high-frequency end. However, the performance of 
the PN approximation deteriorates with decreasing q. 

Effective-One-Body (EOB). The FOB formalism ||5J maps 
the conservative part of the PN dynamics of a compact binary 
system onto the dynamics of a test particle of reduced mass 
ji = mv = OT1OT2/OT in a time-independent and spherically 
symmetric effective metric ds^ff = — A(r; v) dt^+B{r; v) dr^ + 
r^{d9^ + sin^ 9d(p^), which reduces to the Schwarzschild 
metric of a black hole of mass m in the limit v ^ 0. The ex- 
pansions of the EOB potentials A and D = {AB)^^ in terms of 
the Schwarzschild-like coordinate u — m/rare known through 
3PN order as IS [H A = 1 - 2m + 2 V + f 54 _ 41 ^2^ ^ _L 
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^(m^), andD= 1 + 6vm2 + (52 - 6v) Vm^ + ^(m^). To en- 
force the presence of an EOB innermost stable circular orbit 
(ISCO), Ref. ifTSl suggested replacing A by its Pade approx- 
imant of order (1,3), Ap — (1 +au)/{l + hu + cu^ + du^), 
whose Taylor series coincides with the known 3PN result. 

From the recent analysis of slightly eccentric orbits in the 
EOB formalism (V9], the effective-one-body prediction for the 
periastron advance in the limit of zero eccentricity is given by 



^EOB 



D(m)A(m) ' 



(5) 



where = dAp/du, and A=ApA'p + 2u{A'p)^ — uApA'p van- 
ishes at the EOB ISCO. To obtain the invariant relation 
Keob{x), one needs to compute u given x, which we do here 
numerically (for any given v) from the expression of the EOB 
Hamiltonian restricted to circular orbits, and Hamilton's equa- 
tions of motion |fT9ll . The resulting curves are displayed in red 
in Figs.[T}|3] For ^ = 1 and 2/3, the EOB(3PN) prediction ^ 
is within the numerical error up to mQ.(p ^ 0.022. For all the 
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FIG. 1: The periastron advance K of an equal mass black hole binary, in 
the limit of zero eccentricity, as a function of the orbital frequency Clip of 
the circular motion. The NR results are indicated by the cyan-shaded region. 
The PN and EOB results are valid at 3PN order. The lower panel shows the 
relative difference SK/K = (A" — S'nr)/A'nr. 



Other mass ratios, the EOB(3PN) result is within the numerical 
error at all frequencies. When using the EOB potential A(m) 
with 4PN and 5PN terms calibrated to a set of highly accu- 
rate unequal mass non-spinning binary black hole simulations 
II20I . the EOB prediction is within the numerical error at all 
frequencies and for all mass ratios considered. This remark- 
able agreement could be attributed in part to the "pole-like" 
structure at the EOB ISCO in Eq. (|5]l, which is absent from 
the standard PN result Q. 

Perturbation Theory and the Gravitational Self-Force. Ex- 
treme mass ratio inspirals (EMRIs) of compact objects into 
massive black holes, for which m2 ^ m\, are important 
sources of low-frequency gravitational radiation for future 
space-based detectors. Modelling the dynamics of these sys- 
tems requires going beyond the geodesic approximation, by 
taking into account the back-reaction effect due to the inter- 
action of the small object with its own gravitational perturba- 
tion. This "gravitational self-force" (GSF) effect has recently 
been computed for generic (bound) geodesic orbits around a 
Schwarzschild black hole ll2Tl - |23l . In particular, the G'iq) cor- 
rection to the test-mass result /Tschw has been derived f24). 
This calculation determined (numerically) the term p{x) in 



the function W = 1 /K^ = l-6x + qp{x) 



). The re- 



sults are well fitted (at the 10^^ level) by the rational func- 
tion p 14x^(1 + ax)/{\ + )3x + yx^), with a = 12.9906, 
j3 = 4.57724, and 7= —10.3124. (This model improves upon 
the model of Ref. 1 24 1; it is based on a much denser sample of 
GSF data points in the relevant frequency range.) In terms of 
the quantity K we have 
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We used this expression, with the above analytic fit for p(x), 
to produce the dashed blue curves in Figs.[T]-l3] 



FIG. 2: Same as in Fig. [T] but for a mass ratio 9=1 /8. Note that for an 
orbital frequency mQ.ip ~ 0.03, corresponding to a separation r ~ 10m, the 
periastron advance reaches half an orbit per radial period. 



Since p{x) > for all stable circular orbits, the ^{q) GSF 
decreases the rate of precession. Note that the formal diver- 
gence of Kq^p at the ISCO limit (x — ^ 1 / 6) is simply a conse- 
quence of the fact that Q.r vanishes there (by definition), while 
Clip remains finite. This divergence might explain why the 
convergence of the standard PN series seems to deteriorate 
with decreasing q 1.25] , as also illustrated by our results (cf. 
Fig. |3]l. We remind the reader that Eq. (|6]l captures only the 
conservative effect of the GSF, and has a limited physical rel- 
evance near the ISCO, where the actual dynamics transitions 
from an adiabatic quasi-circular inspiral (driven by the dissi- 
pative piece of the GSF) to a direct plunge Il26ll27l . 

We now turn to discuss one of the most striking findings 
of our study. Since q and v = q/{l + q)^ coincide at leading 
order, namely ^ = v + i?{v^), we may recast Eq. ^ as 
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which, unlike Kq^^, is symmetric under mi ^ m2- The solid 
blue curves in Figs. fT||3| show K^^p. Remarkably, while the 
agreement between K^^p and A'nr becomes manifest only at 
sufficiently small q (as expected), Kq^p appears to agree ex- 
tremely well with /Tnr at all mass ratios. This suggests that 
the substitution q v amounts to an efficient "resummation" 
of the <7-expansion, to the effect that much of the functional 
form K{x) is captured by the ^(v) term, even for large q. 

A few heuristic explanations for this behavior may be sug- 
gested, (i) As mentioned earlier, quadratic corrections in v 
enter the PN expression for K only at 3PN [recall Eq. Q], 
and account for less than 1% of at this order. This implies 
that the linear-in-v approximation must be very accurate, at 
least at small frequencies, (ii) The true function K{x;mi,m2) 
must be invariant under exchange mi o m2. The expansion 
in V, Kq^p, satisfies this symmetry by definition of v, whereas 
the expansion in q, Kq^p, does not. (iii) Assuming the co- 
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FIG. 3: The relative difference dK/K = (A" — A'nr) /Ajmr as a function of tlie 
mass ratio q, for mClip = 0.022. Tlie PN and EOB results are valid at 3PN 
order. The shaded area marks the error margin of the NR data. The results are 
qualitatively identical and quantitatively similar for other values of Q.if. 

efficients a„ in the formal expansion K — £„fl„(ji;)v" do not 
increase with n (which, however, only a future calculation of 
higher-order GSF terms could confirm), this expansion will 
exhibit a fast convergence since < V ^ 1/4; the same can- 
not be said of the ^-expansion. 

Comparison of the GSF curves in Figs. [TJ|3] with the NR 
benchmark leads us to yet another important observation. It is 
evident that the second-order GSF correction to K (i.e. the un- 
known term °^ q^) has an opposite sign with respect to the first- 
order term; namely, the second-order GSF acts to increase the 
rate of periastron advance. This is a new result, which illus- 
trates the potential merit of cross-cultural comparisons of the 
kind advocated in this work. 

Summary and Discussion. The advent of precision-NR 
technology allows us, for the first time, to extract accurate in- 
formation about the local dynamics in binary black hole inspi- 
rals (previous studies focused primarily on asymptotic wave- 
forms), and carry out meaningful comparisons with the results 
of analytic approaches to the problem. These comparisons and 
cross-check validations among analytic approximants and NR 
results are crucial for developing faithful analytic waveforms 
to be used in LIGO/Virgo searches. 

Here we focused on a particular aspect of the dynamics, 
namely the relativistic periastron advance. We worked in a 
highly relativistic regime, where the periastron advance can 
reach values as high as half an orbit per radial period (far 
greater than the meagre ^ 43" per century advance of Mer- 
cury's perihelion!) We employed the invariant relation K{Q.,p) 
as a reference for comparison, which is meaningful only in 
the adiabatic regime where the dissipative evolution is "slow". 
For the range of inspiral orbits covered by our NR simulations, 
a measure of adiabaticity is provided by 0.3% < tl^/Q.'^ < 
1.7%. This suggests that inclusion of dissipative effects in the 
PN/EOB/GSF results would not substantially affect our con- 
clusions. The very good agreement between the analytical and 
NR results at low frequency, where the error in /Tnr is small- 
est, also supports this expectation. 



Our direct comparison between perturbative and full NR re- 
sults is the first of its kind. The ff{q) GSF prediction agrees 
with the NR data for small mass ratios (e.g. = 1/8 or 1 /6) 
to within a relative difference of magnitude ^ q^, as expected. 
This provides an extremely strong validity test for both NR 
and GSF calculations. Furthermore, the sign and magnitude 
of the difference /Tnr — ^gsf 8^^^ valuable, hitherto inac- 
cessible information about the second-order GSF effect. 

The above validation test is further reinforced by the 3PN 
result, which shows a good agreement with the NR data at 
small frequencies, or "large" separations (down to r ^ 10m), 
especially for comparable masses (e.g. for ^ = 1 or 2/3). Our 
comparison also reaffirms the expectation that the PN approx- 
imation performs less well in the small mass-ratio regime. 

We find that the EOB(3PN) prediction of the periastron ad- 
vance is in very good agreement with the numerical one across 
the entire range of mass ratios and frequencies considered. 
This result supports the idea that the EOB formalism can de- 
scribe the binary dynamics at all mass ratios. 

Finally, we observe that the simple replacement <7 — > v can 
extend the validity of the GSF approximation far beyond the 
EMRI regime. Indeed, our model Kq^^ agrees very well with 
the NR data at all frequencies, and for all mass ratios con- 
sidered, including the equal mass case. This surprising result 
suggests that GSF calculations may very well find applica- 
tion in a broader range of physical problems than originally 
envisaged, including the modelling of intermediate mass ra- 
tio inspirals, a plausible source of gravitational waves for Ad- 
vanced LIGO/Virgo |28|. 
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